
def mag(counts, zeropoint=32.3):
    """ Calculate magnitude.  Zeropoint is nagnitude of an object with
    1 count per second. """
    return -2.5*np.log10(counts) +zeropoint

fh = pyfits.open('gal.fits')
a = fh[0].data

pixscale = fh[0].header['PIXSCAL1']
pix_per_sqarcsec = 1 / pixscale**2

# find background
#bg = np.median(a,None)
#a0 = a - bg
# tweak the bg after looking at histogram - histogram peak should be
# roughly at 0.
#pl.hist(a0.ravel(),bins=1000,histtype='step')
bg = np.median(a,None)
bg -= 10
a0 = a - bg
#pl.hist(a0.ravel(),bins=1000,histtype='step')
pl.axis('equal')
pl.pcolormesh(a0,cmap=pl.cm.hot)
pl.colorbar()
#pl.pcolormesh(a0,cmap=pl.cm.Paired)

# magnitude of the entire object
mag(a0.sum())

pl.axis('equal')
pl.contourf(a0,10,cmap=pl.cm.hot)
pl.colorbar()
# magnitudes per arcsec**2
mag(pix_per_sqarcsec*800)
mag(pix_per_sqarcsec*1600)
mag(pix_per_sqarcsec*2400)
